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Three dimensional computations of self consistent three species gyrofluid turbulence are carried 
out for tokamak edge conditions. Profiles as well as disturbances in dependent variables are followed, 
, running the dynamical system to transport equilibrium. The third species density shows a signif- 

• icant correlation with that of the electrons, regardless of initial conditions and drive mechanisms. 

' For decaying systems the densities evolve toward each other. Companion tests with a simple two 

fNj , dimensional drift wave model show this persists even if the third species is a passively advected test 

field. Similarity in the transport character of electrons and the trace species does not imply that 
the electrons themselves have a test particle transport character. 



X 



PACS numbers: 52.25.Fi, 52.65.Tt, 52.35.Ra, 52.30.Ex 



I. INTRODUCTION — TRACE SPECIES TRANSPORT 



c/3 ' It is an interesting question whether the electrons and ions in a magnetically confined plasma are transported 
^ . in the same way test particles in a randomly turbulent velocity field would be 0, 0, 0, 0| . One might consider to 
rS i' investigate this by placing test particles or a trace species into an otherwise self consistent turbulence computation or 
experimentally by measuring trace species transport and comparing to the transport of the electrons jBI . Not only the 
O ■ transport coefficients are interesting, but also the basic character, namely, whether or not the transport is a random 
' ^ walk process 0, S H • It might be tempting to conclude that if the transport of the trace species were the same as 
J>-j. that of the electrons, then the electrons themselves should be transported essentially randomly. 

■ It is useful to provide a falsification constraint on this logic. Gradient driven transport in a magnetised plasma 
is usually very different from test particle transport in a fluid, due to the strong coupling between the electrostatic 
potential, which serves as the stream function for the ExB fluid velocity, and the electron density [lol ITll IT^ . The 
coupling is effected by the action of parallel currents and produces a strong correlation between the potential and 

^ ■ density. Ultimately, if the electrons are adiabatic (Boltzmann relation vis-a-vis the potential) then there is no net 
\ transport of electrons even though there is turbulence, which could also be driven by the ion temperature p^ ll4l[T5j . 
i Since ExB transport is ambipolar, it follows that there is also no net transport of ions (neglecting the small charge 
' density implied by a finite ExB vorticity) . Even in robustly electromagnetic tokamak edge turbulence, where the free 
, energy comes from the extent to which the electrons are nonadiabatic, the adiabatic response is always of qualitative 

importance, and there is a strong correlation between the potential and the electron density |l6lll7| . 
, What we will show herein is that in most reasonable circumstances the test particle species, even when it is a tracer 

■ field with vanishing influence upon the electrostatic potential, follows the same dynamics as the electron density. The 
O ' process has a very similar appearance and results ultimately from the same mechanism as "dynamical alignment" 

' between a passively advected density and the vorticity in neutral fluid turbulence |l8j |. It is a manifestation of (1) 
J>-j' the advective nonlinearity being the most effective process in the equation for each quantity, and (2) each of these 
; quantities having a strong direct cascade in its squared amplitude (playing the role of an entropy contribution) despite 
the fact that the flow energy has an inverse cascade 19]. The basic nonlinear dynamics is sufficiently general that these 
j> \ cascade tendencies persist in the presence of strong, dissipative forcing The result is that any initialised imprint 
on the spatial morphology of either variable is quickly transferred to small, dissipative scales, and the subsequent 
. ^ evolution, being the same for both quantities, produces similar morphology especially at scales within the inertial 
range of the turbulence. 

■ " ' The computations are done in three dimensional flux tube geometry, with a fully self consistent three species 
gyrofluid model [U with two ion species plus electrons, generalised to allow nonadiabatic, electromagnetic electron 
dynamics , but restricted to two moments pe r species in the "GEMS" model with clean energetics and clear 
correspondence to fluid edge turbulence models [23 ■ Not only the evolution of the small scale eddies are followed, but 
also that of the zonal profiles (flux surface averages, cf. Ref. Q)- This dynamics is followed to transport equilibrium 
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for cases driven by a source, and for at least one transport decay time for cases with no sources. The trace ion 
species is initialised either the same as for the electrons (random bath disturbances plus a profile), or simply with 
a narrowly localised profile. In some cases both ion species are given equal background densities. In each case, the 
initial differences are eliminated by the evolution of the turbulence, and at late times the morphology of all species 
densities are closely similar, with a high degree of cross correlation between any two of them. 

The main result is that even though the electrons are not transported passively, the density of a trace species is 
transported the same way as the electrons are. A companion test is given by a simple two dimensional drift wave 
model, following only the potential and electron density and their dissipative coupling 11], plus an additional equation 
for a third field which is the density of a strict test particle species. Here, the regime of the electrons response is varied 
from deeply hydrodynamic to adiabatic [TsL [T9l | . In either case, the test particle density closely follows the electron 
density, in all regimes, regardless of whether the electron and test particle densities are similarly or differently driven. 
This should serve as a caution against prematurely concluding that the density in a magnetised plasma are passively 
advected in the case the electrons and a trace impurity species should be observed to display the same transport 
properties. 

Following sections document the GEMS model as used, the results in decaying and driven cases, respectively, and 
the results from the dissipatively coupled, two dimensional model. A discussion and summary is given at the end. 



II. THE THREE SPECIES GEM3 MODEL 



For this study we use the GEMS model, which allows a finite background temperature for all species but does 
not follow the dynamics of the temperatures or parallel heat fiuxes ,2SJ . It allows drift wave and ideal and resistive 
ballooning dynamics according to the parameters 

o 47rpe ( qR\^ ^ Ve TUe ( qR\'^ 



controlling the parallel electron response and 



controlling the perpendicular drift dynamics and the magnetic curvature effects. Both interchange and geodesic 
curvature are retained |24| . Normalisation is in terms of a profile scale L± and the sound speed Cs , giving a natural 
drift frequency of Cs/ L^^. The parallel dynamics is normalised in terms of gi?, where the parallel connection length 
is 2nqR, giving the qR/L± scale ratios in the parameters. The range of scales under consideration is everything 
between the drift scale ps and the profile scale L±. The collisional drift wave regime is roughly C > 1, becoming 
electromagnetic if /3 > 1. Ideal and resistive ballooning enter if f3u}B > 1 or Cujb > 1, respectively. The first 
two conditions are usually satisfied but the latter two are not. Hence, the principal physical process is drift Alfven 
turbulence [l^ [l7| . Geodesic curvature is always important as its role is to regulate the attendant zonal ExB flows 

a. 

The profiles as well as the fluctuations in the dependent variables are followed, so that transport as well as turbulence 
can be computed, even though the equations are still homogeneous in the normalised parameters. Each species 
is characterised by a background charge density, temperature/charge ratio, and mass/charge ratio, given by the 
normalised parameters 

= n^^Z/ue = TjZTe p.^ = Mz/ZMd (S) 

in terms of a reference electron temperature Te and deuterium mass Mo- The scale ratio for the parallel dynamics 
comes in as Cz — fizi^R/ L±)^ . For electrons, ae = = —1 and fie — —(Tte/M^)). For a single component deuterium 
plasma we would have ai — fii — 1 and the nominal ion/electron temperature ratio Ti/T^. With two ion species 
(the second labelled 't') we specify the three parameters independently for each species but then restrict to a neutral 
plasma by taking = 1 — a^. For simplicity we will assume a deuterium/tritium plasma with equal temperatures for 
all species. 

The GEMS model is given by moment equations for the density and parallel velocity of each species. 



dt 



^-BS/i\^+lc(^^G+Tznz) (4) 



-dAll dzUz\\ ~ /~ _ \ / _ ^ 

f^-Qf + ^ ^"^ll ^ yt'C: + TzUz j +JC [ezTzU:z\\) (5) 
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including the effects of the toroidal magnetic drifts also in the parallel velocity moment equation |2lL |25j . The ExB 
advective and parallel derivatives are given by 



4 ^ a_ 

dt dt 



in terms of a Poisson bracket defined by 



dx dy dy dx 



(6) 



(7) 



noting the species-dependent nature of (j)Q. The curvature operator and perpendicular Laplacian are given by 



K- — Slub 



d 



d 



(cos s + sin s) — h sin s — 



dyk 



dx ^ dyk 



dyl 



(8) 



where g^}' — s[s — Sk\ is the off diagonal metric element and s is the shear, and yk denotes the member of the family 
of field aligned coordinate systems which is orthogonal at s = Sk- The gyroaverage operator is defined by 
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(9) 



with pz the species gyroradius in terms of the nominal ps and with k\ the Fourier transform of — (note k\ 
contains the factor The electrostatic potential is determined by the gyrofluid analog 26] of the gyrokinctic 

Poisson equation p^ . 



-pl/2~ , To - 1- 

rg nz ^ < 



(10) 



with Pq always taking the argument (}r^p^. The magnetic potential is determined by the gyrofluid Ampere's law, 



Vi ^11 = Jii 



azUz\\ 



(11) 



hence by the various velocity moments making up the parallel current. 

All of this is in local flux tube geometry defined separately with respect to each drift plane {xyk , defined globally, 
orthogonal at s = s^), each on its own globally field aligned coordinate system, so that the linear parallel derivative 
(d/ds) incurs shifts in the y-direction to reflect the magnetic shear. The simplest version of the geometry is used, 
in which the metric of the coordinate system referenced to each drift plane is unit diagonal on that plane, and the 
magnetic field strength is _B = 1, giving the above Laplacian and curvature operators. This is the "shifted metric" 
treatment; for further detail see Ref. |2^. Elsewhere, the subscript on y^ is omitted (the form of the Poisson bracket 
is unaffected). 

The equations are solved on a domain L Ls, with Ls = 2tt due to the global consistency constraint. The 

boundary conditions on the dependent variables in the drift plane are Dirichlet (vanishing value) at x — Lx/2 and 
Neumann (vanishing derivative) at a; = —Lx/2, and periodic on —Ly/2 < y < Ly/2, with the x-direction implemented 
in terms of the quarter-wave FFT so that the correct value of Pq is obtained for each wavenumber. The boundary 
conditions on s are defined by global consistency f23|, on the interval — tt < s < tt, implemented using the shifted 
metric treatment psj . 

The numerical scheme is one which has been used previously in three dimensional drift wave computation |3Clj |. 
Poisson bracket structures in xy are evaluated with a discretisation which preserves their properties p3l| . The linear 
d/ds and /C terms are evaluated using centred differences. The timestep is a third order scheme using both the 
dependent variables and their time derivatives to build the new dependent variables at the next time step "s^. In 
contrast to most other methods, this scheme is stable to waves. Nevertheless, the direct thermal free energy cascade 
[T^ITgl l creates the need for numerical dissipation at small scale This is effected by a diffusion —v^^d'^/ds^ and a 
hyperdifFusion v±^\ added to each ExB derivative (applying the dissipation to the moment variables but not directly 
to the fields). The spatial grid node count is 64 x 256 x 16 in xys, on a domain given by Ly — AL^ — SOirps — 4L^, or 



simply 4 in normalised units, and 
and runs are taken to approximately t 



27r. The minimum value of kyPs is therefore 0.025. The timestep is r = 0.05, 
2000. The artificial dissipation coefhcients are v\\ = 0.003 and v±^ ~ 0.002. 
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III. GEMS RESULTS FOR DECAYING CASES 



The first set of results under consideration is tlie decaying scenario: the domain is loaded with a density profile and 
density fluctuations, and the system is allowed to relax. Profile maintenance is restricted to a sink term on the right 
hand side of each density equation, 



dt 



dy 



nzp{x) 



p{x) — 0.3 exp 



(12) 



where the integration denotes an average over y. Net transport in the direction \7x into this sink then causes slow 
decay of the overall density profile (zonal average of ne)- 
The initial profiles are given by 



1 



1 — sin — 

Lx 



nto = exp 



{0.1 L^y 



(13) 



with 



Tieo = (1 - at)nio + atuto (14) 

with subscripts {e,i,t} referring to electrons, main ions, and trace ions, respectively, with at the background charge 
density parameter of the trace ions. The electrons and main ions start with a random bath of fluctuations at nonlinear 
amplitude, n, normalised such that the RMS level is 3.0 x S, described elsewhere 113: with 

He = (1 — at)n Ui = n nt — (15) 

The trace ions start solely with their profile. The initial electrostatic potential is then given by Eq. H1U|) . 

The initial evolution is the excitation of "Pfirsch-Schliiter" profile modes, especially the global Alfven oscillation. 
This is caused by the fact that the initial equilibrium is only one-dimensional while the actual one is two-dimensional, 

including Pfirsch-Schliiter currents and flows. It involves the "sideband modes" given by ^J|| coss^ and ^0sins^, 

where the angle brackets denote the zonal average (over y and s). These are excited by the action of the magnetic 
compression K, upon the pressure profile (p), where p = azTzUz, due to the sins term in fC acting upon the ky = 
component; in other words, the geodesic curvature. The decay of the global Alfven oscillation proceeds through 

resistive dissipation, C^Jucoss^ When these sideband modes reach dissipative equilibrium, the Pfirsch-Schliiter 

current is established. Over about the same time interval, < t < 200, the basic nonlinear drift wave mode structure 
is set up ,16i„ .17„ .23i]. Parallel sound wave dynamics, the geodesic acoustic sideband dynamics, and the zonal flows 
all reach statistical equilibrium by about t = 1500 [2J|. For typical edge parameters the profile decay (i.e., transport) 
time scale is comparable to this, so we are dealing with a slowly transient case. However, it does provide a useful 
control against the effects of profile maintenance (driving of the densities using a source) upon the turbulence. 

The standard case, reflecting typical tokamak edge conditions, is given by the parameters S — 0.0159, /3 — 1, 
C — 2.5, Ti — 1, UJB = 0.05, and {qR/L^_Y = 18350. This very roughly corresponds to a physical parameter set of 
Tie = 3 X 10^^ cm^^, Te = 100 eV, i? = 2.5T, g = 3.3, R = 160 cm, and = 4 cm. As noted above, we also use 
deuterium as the main ion species and tritium as the trace species, with Tt — I and /it — 1.5 and Mu/nie = 3670, 
and either a* = (strict test particles) or at = 0.1 (minority species). 

Selected time traces of the case with at = are shown in Fig.^ In the left hand frame, the signal An is the grid node 
average of azTzn1/2, giving the total thermal free energy, and Ap is 0^/2, tracking mostly the zonal flows, and 
A^, is ri^/2, tracking mostly the activity of the turbulence, and Ff, is rieV^, the transport, with = —S{d4)/dy). In 
the right hand frame, the signal is az4>G^z/2, the ExB energy, and Eb is azA\]Uz\\l'2,^ the magnetic energy. 

Here, Vl — —^^azriz defines the generalised vorticity (at zero FLR it is just the ExB vorticity). The signature of 
Ef. and Eb is that of the global Alfven oscillation. The equilibrium is well established after about t — 250. The 
other traces show that the overall statistical equilibration time of the turbulence and zonal fiows is comparable to the 
transport decay time, and that the vorticity and transport track each other closely. By t = 1500 the total amount of 
electrons, n,, integrated over the volume, has decayed by about half. 

The evolution of the decaying profiles (zonal averages, over all y and s for particular x) of the electron and trace 
densities is shown in Figs.[21and|3| In these units, the maximal initial value is unity, and the first frame at i = 10 
shows that the turbulence is already deconstructing the initial profile of nt. The trace species distribution spreads 
until the value at a; = —Lx/2 (—32 in units of ps) rises from zero at t = 100. Thereafter the trace profile fills in as 
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the electron profile decays. When the value at x = reaches that at a; = —Lx/2, the nt profile itself begins to decay. 
At late times, the profiles of and are similar. As the driven cases discussed in the next Section will show, this 
behaviour is not indicative of a particle pinch in either species. 

The most prominent amplitude/energy spectra (squared amplitudes, Fourier decomposed y ^ ky, averaged over 
< X < Lx/2 and all s) are shown in Fig. 0] averaged over the intervals 500 < t < 600 and 1500 < t < 1600, 
respectively representing the early and late stages of the turbulence. These are reflective of typical drift wave mode 
structure as we expect in this regime; note especially that as a function of ky (labelled n) follows cjP (labelled p) 
throughout the spectrum, as the average value of the parallel wavenumber fc|| self adjusts to whatever is necessary 
for the parallel currents to balance the nonlinear excitation rates of the ExB vorticity. The vorticity spectrum itself, 
fl^ (labelled w), is much flatter, extending all the way to ky = 1 in units of p~^. The transport spectrum, Uev'^, also 
shown for both early and late stages, peaks at relatively low ky but is broadband. We flnd similarly standard drift 
wave mode structure results for the parallel envelopes and transport and energy transfer spectra, as presented in Refs. 
[l3,|23, not shown here. 

The spatial morphology at t = 50 and 400 is shown in Fig. [SJ as contour patterns in the xy-plane at s = 0. The 
form of nt and rie start out very differently owing to the initialisation. The spatial morphology however becomes 
quickly similar after the turbulence deconstructs the original proflle of nt (transfer free energy from fc^ = to the 
entire /cy-spectrum). At i = 50 the part of the contour distribution for a; > is already very similar, even though the 
values at a; = are still close to zero. However, the spectra and even the spatial morphology evolve toward each 

other on all scales, including the profile scale. After t — 400 they remain very close. The fact that the spectrum of nt 
is slightly flatter on the smallest scales is reflected in the somewhat noisier form of the morphology of nt , seen most 
clearly in the outermost contour. 

The closeness of nt and Tie may be quantifled by their cross correlation and phase shift distributions. These are 
more usually em ploy ed to evaluate the relation between ng and 0, to distinguish between drift and MHD mode 
structure [Ta. [l7ll23| ). An "adiabatic" relationship is characterised by a close cross correlation and a narrow phase 
shift distribution peaked near zero. A "hydrodynamic" relationship, the one expected of passive scalar advection in 
a neutral fluid, is characterised by a very weak or zero cross correlation and a wide phase shift distribution centred 
upon 7r/2, the value which gives energetically maximum gradient drive. The fields in question are sampled on the 
xy-plane at s = 0, with the ky = mode stripped out, over the interval 500 < t < 1500 with the turbulence saturated 
and still well developed, and for 5 < a; < 27 to keep clear of the outer boundary sink region and remain in the 
gradient region for nt- These diagnostics are shown for both {rie,0} and {nt,ne} in Fig. The first pair show 
close correlation (overall: 0.721) and low phase shifts (positive, narrowly distributed, below 7r/4) indicative of drift 
wave mode structure. The second pair shows extremely tight correlation (overall: 0.937) and a narrow phase shift 
distribution centred upon zero. 

The expected form for a passive scalar quantity completely uncoupled to the potential would be hydrodynamic, 
with random phase shifts in a case with no gradient. By contrast, although in the dynamical equations there is no 
direct coupling between ng and nt, they are very closely correlated and their phase shift distribution is narrow and 
peaked at zero for every nonzero ky in the spectrum. Taking the two combinations from this figure together, we find 
that the electron and trace ion density fluctuations track each other very closely, despite the fact that the electrons 
themselves are not in a passive advection relationship to the ExB eddies. This is the phenomenon we refer to as 
"dynamical alignment" and the center of the phase shift distribution (essentially zero) and the value of the overall 
cross correlation (larger than 0.9) give it a quantitative basis. 

Testing against the presence and absence of a back reaction of (j) to nt, this basic case with (3=1 and at — 
was run again with at — {0.1,0.2,0.3,0.5}. The same results as presented above for at = were found, with the 
cross correlation values for at — 0.1 given by 0.763 for n,. and (j) and 0.947 for nt and Hg, respectively. Cases with 
/? = 5 with at = and 0.1 were also run, looking for changes closer to the MHD regime in which the mode structure 
should be more hydrodynamic. The basic mode structure shows the changes resulting from the stronger role of the 
vortices at kyPs — 0.05 to 0.1 (convcctive cell modes on the scale of L±^, which in the fc^-spectra sit higher above the 
broadband drift wave turbulence and show much larger phase shifts. Nevertheless, the same closeness of nt to rie is 
found. The cross correlation and phase shift information as in the nominal case are shown in Fig. [7| for the case with 
/3 = 5 and at = 0.1. The fact that this case is in the MHD regime is shown by the hydrodynamic relationship between 
rie and 4'] the overall cross correlation value is —0.0259. The persistence of the dynamical alignment is shown by the 
relationship between nt and ng, which gives zero phase shifts and a overall cross correlation of 0.986. 

Regardless of the closeness of the relationship between ng and 0, which varies with parameters (in this case (3 and 
C), the correlation between Ug and nt is close to unity and the phase shifts are close to zero. Clearly, nt is following 
ng regardless of the latter's relationship to (/). 

What we now have to do is to determine the extent to which this result might be a fortuitous consequence of 
how the model is set up. To this end, we compare to driven cases run to transport equilibrium in the next Section, 
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and then evaluate the phenomenon more fundamentally by studying the simplest dissipative drift wave model in the 
adiabatic and hydrodynamic limits, in the Section after that. 



IV. GEM3 RESULTS FOR DRIVEN CASES 



We now turn to cases in which the profiles are maintained by sources so that the runs may be carried for arbitrary 
time. The setup is exactly as for the decaying cases, with the addition of the drive term 

^ = --- + S.{x) (16) 

for each density variable. Typically we use Si — Se so as not to drive vorticity on the inner boundary. The basic 
profile source term has the shape of a Gaussian with 1/e- width O.^L^ centred upon x = —Lx/2 and a time constant 
chosen to anticipate the transport time scale. Narrower main source profiles (width O.lLx) with correspondingly 
larger coefficients (i.e., at constant total source) were found to excessively excite artificial MHD modes. The trace ion 
source was centred upon x = to set up a gradient region for a: > and to look for a pinch (finite gradient with no 
interior source) in the x < region, so rito was used as a template. The source terms were chosen as 



Se^ Si ^ Se exp 



St = stnto (17) 



(0.3L,)2 



for electrons, main ions, and trace ions, respectively. The parameters Se and St are free; nominally they are both set 
to 10-3. 

The standard case is the same as the one from the previous Section, with = 0, so that the trace ions are proper 
test particles. The same time traces are shown in Fig. |S1 which can be compared to Fig. ^ The amplitudes and 
transport, in the left hand frame, show that the dynamics is well saturated after t = 1000, so that late time averaged 
diagnostics are taken over the interval 1000 < t < 2000. The time traces of the ExB and magnetic energies, labelled 
Ee and Eb, respectively, are shown in the right hand frame, indicating that the global Alfven oscillation damps away 
well before the saturated stage. The short time scale visible in the Ap signal is indicative of the geodesic acoustic 
oscillation. The long time scale oscillation from the decaying case is absent in this figure, showing that it is a response 
to the overall profile decay in the other case. 

The evolution of the profiles of the electron and trace densities in this driven case is shown in Figs. El and 1101 
measured the same way as in Figs. |21 and |21 for the decaying case. In these units, the maximal initial value is unity, 
and the first frame at t = 10 shows that the turbulence is already deconstructing the initial profile of n*. Up until 
t = 100 the evolution is the same as in the decaying case. The trace species distribution spreads until the value at the 
left boundary {x = —32 in units of ps) is as large as that at the source {x = 0). Transport equilibration, the temporal 
convergence of these profiles, sets in after t = 1000. At late times, the profiles of and rit are similar outside the 
source regions. The flatness of nt for a; < is indicative of the absence of an impurity pinch. For this we note however 
that temperature disturbance and gradient dynamics is not followed, so there remains the possibility of an impurity 
pinch driven by the ion temperature gradient. Nevertheless, the background temperatures of the three species are 
equal, so this does show that the existence of warm ions does not by itself lead to an impurity pinch nor does the 
trace and main ion temperature lead to a "neoclassical" pinch, the underlying dynamics of which would indeed be 
contained in this model (geodesic curvature acting on the gradients through diamagnetic fluxes, and the interaction 
of this with parallel flows in the axisymmetric part of the dynamics). 

The spectra of the amplitudes and transport are shown in Fig. 1111 averaged over the saturated stage, given by 
the interval 1000 < t < 2000. These show the amplitudes of He and (p peaking at large scales {ky < 0.1 in units of 
pj^), while the transport is broadband almost all the way to fcj, = 1. The vorticity is flat to fcj^ = 1. All of these are 
indicative of drift wave mode structure. The spectrum of rit follows Hg. 

The spatial morphology at i = 2000 is shown in Fig.^l as contour patterns in the xy-plane at s = 0. The form of 
rie and rit are closely similar. The fact that the spectrum of nt is slightly flatter on the smallest scales is reflected in 
the somewhat noisier form of the morphology of n(, seen most clearly in the outermost contour. From time to time a 
small scale structure appears in nt which is not reflected in n^, but the occurrences are rare enough not to affect the 
statistical result. _ 

The cross coherence and phase shifts between n^ and (j) and between rit and rie are shown in Fig. 1131 sampled over 
the saturated stage the same way as in the decaying case. The first pair show that the signature of drift wave mode 
structure (overall cross correlation: 0.674) is not affected by the proximity of the source. For realistic drive levels 
there is therefore no qualitative difference between "ambient" and "source driven" turbulence, provided the physical 
separation of scales is well reproduced in the computations (note that a transport equilibration time of 1000 with 
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tokamak edge values of L±/cs is already very short). The second pair show that the presence of the source just to the 
left of the measured region does not affect the dynamical alignment result (overall cross correlation: 0.920), as can be 
judged by comparing with Fig. El But if this is measured over the whole radial domain, then the effect of not only 
the source but also the region interior to it is found to be large. So the conclusion of dynamical alignment between 
these two disturbance fields does depend on the relative strength of the nonlinear dynamics. 

Having now determined that the three dimensional model recovers the same basic result in source free gradient 
regions, whether or not it is driven or allowed to decay, we now turn to the simplest two dimensional dissipative drift 
wave model for further comparison. 

V. RESULTS IN A SIMPLE 2D DRIFT WAVE MODEL 

We distill the central physics behind the above results by comparing them to those from a simple two dimensional 
Hasegawa Wakatani drift wave model 0, augmented by a continuity equation for a test species. The equations are 

^+^E-V^^D{^~n^-lC{ne) + Vd4> (18) 
^+VB-V(ne + ne)-i:'(0-?ie) + K, - (19) 

H- + V£ • V (/ + /) = (20) 

for disturbances in the ExB vorticity, electron density, and test particle density, respectively. The domain is a single 
xy drift plane, doubly periodic. Background gradient drive terms appear for both the electrons (jig) and test particles 
(/); their x-derivatives yield the background gradients. The ExB advective derivative and vorticity are given by 

v.-v=f|?|:-^^) n = vl~^ (21) 




dx dy 

The dissipative coupling parameter D serves as a model for the adiabatic response. The hydrodynamic and adiabatic 
limits are D ^ and D oo, respectively. For D = the electrons and test particles have the same equation. 
For D > 1 they are obviously very different. A simple interchange curvature model with /C = cuBd/dy and a long- 
wave damping coefficient are included to study saturated interchange turbulence. Note that this model uses the 
traditional "gyro-Bohm" local normalisation, with the factor of S folded into the amplitudes. 

For the purposes of this test series, the zonal flow question is avoided by setting temporal changes to the ky = 
component of ft to zero. The Colella 'ss'l/Van Leer's^ MUSCL scheme is applied to d/dt as in previous versions 
of the 3D electromagnetic gyrofluid model 22] (the linear drive terms are combined into this, e.g., b y d efining N = 
TLe + Tie). The terms involving D are done via the implicit scheme of previous 2D drift wave models [33 • In this, D 
is applied only to the ky ^ components. The background profiles are chosen differently, with due/dx = — 1 and 
df /dx — ~LOz cosTTx/Lx so that Tie and / are driven differently regardless of the value of the gradient parameter w^. 
The domain size is given by L^, — Ly = 2'it/K, with K nominally set to 0.1 along with nominal values of D = 0.1 and 
LOz = 1. There are 64 x 64 grid nodes for all cases. The timestep is set to r = 0.05. All runs which reached saturation 
were taken to t — 1000 with the interval 400 <t< 1000 as the saturated stage. 

The variations for the drift wave cases with wb = Vd = ^ were 
D = {0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1., 3.} with nominal K and w^, K = {0.025, 0.05, 0.1} at = 0.01 and nominal 
Wz, and u)z = {0.1,0.2,0.5, 1.} at nominal D and K. The cases with D = 0.001 and 0.003 did not saturate, as thin 
radial flow jets formed and simply grew, so another scan with D — {0.001,0.003,0.01} was done aX K ^ 0.025 and 
nominal lUz- In all cases in which saturation was found, the cross correlation between (j) and Ue was larger than 0.85. 
This is due to the nature of the turbulence in this model: at low D the spectrum migrates to large enough scale that 
D can compete with k\uj, with uj close to the diamagnetic value of ky, thereby guaranteeing strong adiabatic coupling 
for the long-wave side of the energy containing range. At high Z), with all energetic transfer phenomena becoming 
weak at large scale, the turbulence again migrated there, and with the weak nonlinearities the cases for D > 1 are 
strongly adiabatic with most of the energy again at long wavelength; the D = 3.0 case was taken to t = 4000 with 
the saturation stage taken as the 2000 <t < 4000 interval. All of this has been the subject of previous study [llll20| . 
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Regardless of any of this behaviour, the main point of this present study remained. The cross correlation between 
Uf, and / was measured at 0.984 or larger for all cases. For D — 0.01 and larger at nominal K and ujz these correlation 
values were 

{0.999,0.998,0.992,0.985,0.984,0.987} respectively. For co^ = 0.1 and larger at nominal D and K all values were 
identical at 0.992. For K = 0.025 and larger at D = 0.01 and nominal they were {0.998, 0.9994, 0.999}. For 
D = 0.001 and larger at X = 0.025 and nominal lu^ they were {0.99995, 0.9998, 0.998}. The cross correlation and 
phase shift distributions between ng and 4> and between / and ng are shown in Fig. 1141 for the nominal case. The 
first pair show the expected near-adiabatic drift wave mode structure and the second pair show very close correlation 
between the densities. 

The purely hydrodynamic limit is found by setting D — and then ujb — 0.1 to drive interchange turbulence and 
= 0.01 to allow it to saturate. The cross correlation and phase shift distributions for this case between rie and (j) 

and between / and ng are shown in Fig. El The first pair show the expected hydrodynamic mode structure and the 

second pair show very close correlation between the densities. 

The phenomena of dynamical alignment therefore remained for all cases, surviving the control check represented 

by this simple 2D drift wave/interchange model. 

VI. CONCLUSIONS 

A set of three dimensional computations of turbulence in an electromagnetic drift wave regime commensurate with 
tokamak edge turbulence of a fully self consistent, two-component plasma has been given. The principal focus has 
been on the second ion species as a proper test species, namely, with no back reaction on the polarisation, as controlled 
by a background charge density parameter nominally set to zero. Variation of this parameter was found to produce 
negligible changes. Both decaying and source driven cases were considered. This difference also produced no difference 
in turbulence regime, though with the sources present the turbulence and transport were able to fully saturate, leaving 
more activity at very small scale at late times. A set of two dimensional computations within a dissipative drift wave 
model with a test species continuity equation was also given as a control case. The numerical schemes and boundary 
conditions in the 2D and 3D models were different. Nevertheless, the same conclusion with regard to the correlation 
of the test ion species was always the same: the test ion and electron densities were found to be correlated to much 
better than 80 percent in all cases, and better than 90 percent in most, as measured outside the source region in the 
3D driven cases. 

It is important to note that the relationship between the electron density itself and the electrostatic potential (the 
stream function for the ExB eddies) was itself quite different among some of these cases. The adiabatic response of the 
potential back to the electrons is parameter dependent, with the electrons usually not in a test particle relationship 
vis-a-vis the eddies. Nevertheless, so long as there is a significant nonadiabatic electron density response, the test 
particle transport was found to conform to whatever the electron transport result was. 

Indeed, these result suggest that the character of the adiabatic response is not the electrons following the potential, 
but the other way around: the potential pushes all species similarly through ExB advection, but the form of its 
disturbances and hence the ExB eddies follows that of the electron pressure disturbances due to the adiabatic response. 
The similarity of transport of all species regardless of fractional composition now becomes an expected result, since 
the ExB advection is the most robust effect in any of the continuity equations. The adiabatic response is therefore 
mostly a matter of the parallel current mediating the dynamics of polarisation. 

There is an important experimental implication in this: Experiments are underway to measure the transport of 
trace tritium ions in tokamaks |5|. Transport of the trace ions is to be compared to that of the electrons. It is very 
important to note that if the two transports are found to have the same character, it does not follow that the electrons 
themselves transport like test particles, namely, that they should follow a random walk process. The existence of 
adiabatic coupling makes the electron transport different from this. The results in this work show that the potential 
finding that trace ion transport would be similar to the electron transport does not by itself give an indication of the 
character of the electron transport itself. 

It remains to establish what these dynamical relationships become in a regime where the passing electrons are 
completely adiabatic, and particle transport is solely due to trapped electrons, and thermal transport is mostly due to 
the ions, through the standard ITG scenario which receives most current attention 15] . Moreover, the question of an 
impurity driven particle pinch remains open, since the actual temperature fiuctuation and gradient dynamics is not 
included in the GEM3 model. Future manifestations of this work will address these questions with a more detailed 
model designed to treat them in either edge or core plasma regimes. 
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FIG. 1: Time traces of free energy amplitudes and transport, showing initial saturation of the turbulence after t — 500 and 
overall decay on a scale of about 1500, all in units of L±/cs. 
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FIG. 2: Profiles of the electron (ne) and trace ion (nt) densities at early times i = 10 and 100 and 300 in the evolution. The 
initial value is unity in these units, and x is in units of ps- The trace ion profile is deconstructed by the turbulence and fills in. 
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FIG. 3: Profiles of the electron (ric) and trace ion (rit) densities at medium to late times t = 500 and 1000 and 2000 in the 
evolution. The initial value is unity in these units, and x is in units of p,. When the value of the trace ion profile at the left 
boundary becomes the maximum, both profiles decay together. 
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FIG. 4: Spectra of selected quantities at early {t = 500) and late (t — 1500) stages of the turbulence, with ky in units of 
pj^. Squared amplitudes of (f> and rie and rit and the vorticity fl are labelled by 'p' and 'n' and 'z' and 'w' respectively. The 
tendency of and rie to follow each other with fl much flatter is a basic signature of drift wave turbulence. The electron ('n') 
and trace ion ('z') transport fluxes also follow each other. 
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FIG. 5: Spatial morphology of the electron (n^) and trace ion (nt) densities at f = 50 and 400, during saturation. The contour 
interval is shown in units oi 5 = 0.0159, and x and y are in units of ps in each case. Starting very different, their spatial forms 
evolve towards each other. 
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FIG. 6: Cross correlation and phase shift distributions between the electron density (rie) and electrostatic potential (0), and 
between the trace ion density (nt) and rie, averaged over the interval 500 < t < 1500, with ky in units of pj^ , measured as 
described in Refs. [iM Il7|. Note that (j) serves as the stream function for the turbulent ExB eddies. The first pair show the 
usual signature of drift wave mode structure (close correlation and small phase shifts at all wavelengths), while the second pair 
show the dynamical alignment between the trace ion and electron densities. This is the case with at = 0; the one with at — 0.1 
was quantitatively very similar in all respects, as detailed in the text. 
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FIG. 7: Cross correlation and phase shift distributions between the electron density (rie) and electrostatic potential (0), and 
between the trace ion density (nt) and n^, for the case with at = 0.1 and ~ 5, averaged over the interval 500 < t < 800, with 
ky in units of pj^ . The first pair now shows the high-beta ideal ballooning mode structure (very low correlation and phase 
shifts near 7r/2). The second pair shows the dynamical alignment between the trace ions and the electrons. Compare to Fig. 
|S| The dynamical alignment persists regardless of the actual character of the electron transport itself. 
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FIG. 8: Time traces of free energy amplitudes and transport, showing initial saturation of the turbulence after t — 500 and 
equilibration of the zonal flows and transport after about t = 1000, all in units of L±/cs. 




FIG. 9: Profiles of the electron (rie) and trace ion (nt) densities at early times t = 10 and 100 and 300 in the evolution. The 
initial value is unity in these units, and x is in units of ps- The trace ion profile is deconstructed by the turbulence and fills in. 
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FIG. 10: Profiles of tfie electron (we) and trace ion (nt) densities at medium to late times t — 500 and 1000 and 2000 in the 
evolution. The initial value is unity in these units, and x is in units of p^. Transport equilibration occurs after about t — 1000. 




FIG. 11: Spectra of amplitudes and transport in the saturated stage 1000 < t < 2000, with ky in units of p7^- Squared 
amplitudes of <j) and rie and nt and the vorticity Q are labelled by 'p' and 'n' and 'z' and 'w' respectively. The tendency of </> 
and Tie to follow each other with much flatter is a basic signature of drift wave turbulence. The electron ('n') and trace ion 
('z') transport fluxes also follow each other. 
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FIG. 12: Spatial morphology of the electron (n^) and trace ion {nt) densities at two different times near the end of the run. 
The contour interval is shown in units of 5 = 0.0159, and x and y are in units of ps in each case. Starting very different, their 
spatial forms track each other closely in saturation. 
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FIG. 13: Cross correlation and phase shift distributions between the electron density (rie) and electrostatic potential (<^), and 
between the trace ion density (nt) and ne, for the driven case with = 0.1, averaged over the saturated stage 1000 < t < 2000, 
with ky in units of pj^, measured as described in Refs. UTOliITJ . Compare with Fig. |S| Provided the sampling is done exterior 
to the source region, the results are the same. In particular, the drift wave mode structure signature remains whether or not 
the turbulence is "ambient" or "source driven" . 



20 



-4 



1 1 1 


, , , 




, , , 


, , , 



-4 
4 



-4 



-4 




— 7T 



a TT 



1 1 1 1 1 1 


K 




1 ' 






- 










P 


1 1 1 1 1 1 1 


10-' 


1 


1 1 



— 7T 







FIG. 14: Control check using the nominal case from the 2D dissipative drift wave model with a test particle density field /. 
Cross correlation and phase shift distributions, with ky in units of pj^, between the electron density (n^) and the electrostatic 

potential (0), and between / and rie. The first pair shows the drift wave mode structure and the second pair shows the 
dynamical test density alignment. 



21 



-4 



— I 1 1 1 1 1 r 




J I I I I I L 



-4 
4 



-4 



1 1 1 


1 1 1 


■ / 

, , , 


, , , 



-4 




— 7T 



a TV 




— 7T 







a TV 



FIG. 15: A further control check using the 2D interchange model with a test particle density field /. Cross correlation and 
phase shift distributions, with ky in units of pj^ , between the electron density (ne) and the electrostatic potential ((/>), and 

between / and rie. The first pair shows hydrodynamic mode structure and the second pair shows the dynamical test density 
alignment. 



